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Ph| abstract 

Q ! Observations of high-redshift Ly-a sources are a major tool for studying the high-redshift 

^ ■ Universe and are one of the most promising ways to constrain the later stages of reioniza- 

^ I tion. The understanding and interpretation of the data is far from straightforward, however. 

. We discuss the effect of the reionizing intergalactic medium on the observability of Ly-a 

sources based on large simulations of early structure formation with radiative transfer. This 
takes into account self-consistently the reionization history, density, velocity and ionization 
^ . structures and nonlinear source clustering. We find that all fields are highly anisotropic and 

■ as a consequence there are very large variations in opacity among the different lines-of- sight. 
\ The velocity effects, from both infall and source peculiar velocity are most important for the 

■ luminous sources, affecting the line profile and depressing the bright end of the luminos- 
CN I ity function. The line profiles are generally asymmetric and the line centers of the luminous 

sources are always absorbed due to the high density of the local IGM. For both luminous and 
average sources the damping wing effects are of similar magnitude and remain significant 
\ until fairly late, when the IGM is ionized between 30% and 70% by mass. 

. The ionizing flux in the ionized patch surrounding a high density peak is generally 

strongly dominated, particularly at late times, by the cluster of faint sources, rather than the 
central massive galaxy. Our results reproduce well the observed mean opacity of the IGM at 
z ^ 6. The IGM absorption does not change appreciably the correlation function of sources at 
' high redshift. Our derived luminosity function assuming constant mass-to-light ratio provides 

5^ I an excellent match to the shape of the observed luminosity function at z = 6.6 with faint- 

end slope of a = —1.5. The resulting mass-to-light ratio implies that the majority of sources 
responsible for reionization are too faint to be observed by the current surveys. 

Key words: cosmology: theory — diffuse radiation — intergalactic medium — large-scale 
structure of universe — galaxies: formation — radio lines: galaxies 
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1 INTRODUCTION 

The reionization of the universe was the last global transition of 
the Intergalactic Medium (IGM), from fully-neutral after cosmic 
recombination at z ^ 1100 to fully-ionized as we see it today, 
caused by the radiation from the first stars. Currently there are still 
only very few direct observational constraints on this epoch. The 
lack of Gunn-Peterson trough in the spectra of high-redshift sources 
indicates a low mean neutral fraction xhi < 10~^ out to redshift 
z ^ 6, which implies overlap was achieved sometime before that, 
at 2;ov > 6 . 

^On the other hand, the WMAP 3 -year data dSpergel & et all 

l2007h yielded a fairly high value for the integrated Thomson elec- 



tron scattering optical depth to the surface of last scattering, at 
Tes = 0.09=b0.03. This requires a significant ionized fraction out to 
high redshifts, z > 12, and thus implies an extended reionization. 
The optical depth by itself does not put very stringent constraints 
on the possible reionization histories, however. The reas on for this 
is the self-regulated nature of the reionization process Clliev et al.l 
l2007ah . whereby the Jeans-mass filtering of low-mass sources in 
the ionized regions results in the Tes and Zov, the overlap redshift, 
being only loosely related. The overlap redshift Zov is determined 
by the abundances and efficiencies of the high-mass sources, whose 
formation is not suppressed by reionization, while Tes depends on 
both high- and low-mass sources. Thus, varying the ionizing effi- 
ciencies of the small sources yields a wide range of Tes values for 
the same value of Zov. 
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dramatically in the near future due to a number of large ob- 
servational projects which are currently under way. The 21- 
cm data from high redshifts contains potentially the richest set 
of information since the signal is inherently three-dimensional , 



on the sky and in red shift/frequency (e.g. *Mada u et al 



iTozzi et al 



200Q|; llliev e t al. 2002, 2003; Zaldarriag a et al 



1997; 



2004 



'Mellema et al." 2006b; Shapiro et all i2006 : llliev et al] [2008 a. 
f urlanetto et al. (2006a) for a detailed recent review). The fea- 
tures that could be derived include the full reionization history, 
geome t ry, statistics and individual bright feature s dMellema et al.l 
l2006bl : iFurlanetto et al.ir2006al : llliev etaDl2008ah . There are sig- 
nificant challenges to be overcome, however, particularly related 
to precise subtractio n of the very strong foregrounds present at 
low frequencies (e.g. IFurlanetto et al.ll2006al : iMorales et al.ll200^ : 
llliev et al.l liooia'). 

The patchiness of reionization creates secondary temperature 
anisotropics i n the CMB thro u gh the kin e tic Sunyaev-Zel'dovich 
effect (.Gruzinov & Hul Il998l: IHuI l2000l: iGnedin & Jaffd l200ll : 



ISantos et al.l2003l : lMcOuinn et al.^2005[: llliev et al.l2007bh. as well 



as polarization anisotropics (Hu 2000": Haiman & Holdej l2003l : 
ISantos et al. 2003; Mortonson & Hu 2007; Dore et al. 200%. Un- 
like the 21 -cm signal, the reionization signatures in the CMB are in- 
tegrated over the reionization history and contain no frequency in- 
formation. However, the typical scales of reionization a re reflected 
in a ch aracteristic peak of the kSZ anisotropy signal ( Ili ev et all 
l2007b h. and the shape of the power spectrum is dependent on the 
reionization parameters (source efficiencies and small-scale gas 
clumping). CMB anisotropy observations can therefore provide us 
with key information about the process of reionization and since 
its systematics are different it would be an important complement 
to the 21 -cm studies. One could also combine these observations 
more directly, by using 21 - cm observations to derive the Thomson 
optical depth fluctuations ('Holde r et all l2007h . Small-scale CMB 
anisotropy and polarization measurements would be quite difficult 
due to the weakness of these signals, but are within th e abilities of 
modern detectors Clliev et al.i2007bl : rDore et 

Narrow-band searches for high-redshift Ly-a emitters 



redshifts, currently up to z ^ 7 (e.g. 


Hu et al.l l200i l2005l: 


Kodaira & et 


al.l l2003l: iTam 


suchi & et al.l 


2005: Rhoads & et alJ 


20031: IStanwav & et al.l 


20041: Malhotra & 


Rhoads" l2004l: 
Kashikawa et alJ 


Bunker et al. 


l2006l: Shimasaku et al.l l2006l: 


20061: lOta et al.ll2008). Together with studies of the Lv-a resonant 



lFa^l2004 iFan & et al.ll2006h they provide important indep endent 
approaches for studying reionization (see iFan et al.l l200d for a 
recent review). The optical depth for hy-a re sonant absorption 
in neutral IGM a t high redshifts is quite large (IShklovskiilll"964l : 
iGunn & PetersonI 1 19651) . thus absorption studies are mostly 
sensitive to very low hydrogen neutral fractions, typically below 
XHi ~ 10~^, otherwise the absorption saturates. This technique 
is thus best suited for studying highly-ionized regions and the 
very end of reionization and has been very successful for setting 
low limits for the redshift at which reionization was completed. 
On the other hand, hy-a emitter surveys do not require a very 
low average hydrogen neutral fraction and t hus in p rinciple can 
probe further into the reionization epoch jMalhotra & Rhoads 
[2004h . Othe r related probes inclu d e damping wing m e asure- 
ments (e.g. iMiralda-Escudd 1 19981 : iMesinger & HaimanI l2004l : 



iTotani et all l2006h. quasar HII region sizes and their evolution 



dark gaps in high-z spectra fe.g. lFan & et aP toodlGaTleranietaLl 
l2008h . 



The correct interpretation of these data for reionization is far 
from straightforward, however. The strong dependence of Ly-a ab- 
sorption on the neutral fraction and gas density means that both 
should be modelled with certain precisi on. High-redshift sources 
are rare and strongly c lustered (see e.g. iBarkana & Loebl l2004bl : 
IFurlanetto et al.ll2004bl : llirev et al . 2007ah . As a result, H II regions 
generally contain multiple ionizing sources and grow much larger 
than the ones created by individual sources. This minimizes the ef- 
fect of the damping wings and increases the transmission, allowing 
the detection of more and fain ter sources than would be naiively 
expected ( Wvithe & Loebl2005b . However, while this is the generic 
expectation, the actual effect would be dependent on the exact ge- 
ometry of reionization, e.g. sources close behind a neutral region 
will be damped even if they are inside a very large H II region. Sim- 
plified models typically assume spherical ionized regions and either 
ignore source clustering or assume linear bias (e.g. Cen & Hai manI 
I2OOO; Santos 2004; Furlanetto et al. 2004a; Haiman & Cen 200^. 
Vyithe et al.. i2005i : iKramer et al.i i2006 : iWyithe & Loeb 2007). In 
practice the large ionized regions form by local percolation of 
multi ple smaller one s and as a consequence are highly non spher- 
ical ( l Uev et al.ll200^ : lMellema et al.ll2006d : llliev et al.l20^ a). The 
inhomogeneous cosmological density fields and non-equilibrium 
chemistry effects (particularly in recently-ionized gas) further com- 
plicate the picture and point to the need of following the cosmolog- 
ical structure formation and the reionization history of a given re- 
gion. A proper account of all these effects can only be done through 
detailed cosmological radiative transfer simulations. 



A number of radiative transfer methods have been developed 
in recent years and now they are reachin g a certain level of m aturity 
and are producing fairly reliable results (llliev & et IDl2006h . How- 
ever, performing large-scale reionization simulations, as required 
for Ly-a studies is still technically very challenging. Recently 
iMcOuinn et al.l (120071) used large scale structure formation numer- 
ical simulations postprocessed with radiation transfer to study the 
observability of Lyman-a emitters at high redshifts and what these 
can tell us about reionization. In order to achieve high dynamic 
range, these authors employed a subgrid model for the collapse of 
the smallest halos. Another, se mi-numerical approach was used by 
iMesinger & Furlanettd (l2008ah . who took the linear density and ve- 
locity fields at early time (essentially the initial conditions for an 
N-body simulation) and used an excursion- set approach combined 
with a first-order Lagrangian theory to "paint" the H II regions on 
the density field. This procedure provides large dynamic range at a 
low cost, but at the expense of making significant approximations 
and thus cannot fully replace full numerical simulations. 



In this paper we use the results of large scale numerical sim- 
ulations to study the transfer of Ly-a through the IGM. In § 2 we 
briefly describe our simulation method. In § 3 we describe the evo- 
lution and environment of a rare, massive source, similar to the ones 
that are currently observed. In § 4 we address the observability of 
Lyman- a emitters, considering the reduction of the transferred line 
flux due to the absorption in the IGM and luminosity functions. In 
§ 5 we describe the effects of the patchiness of reionization on the 
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2 SIMULATIONS 

Our simulation results follow the full, self-consistent reionization 
history in a large volume of (100 h~^Mpc) ^ and were desc ribed in 
detail in IHev et al. (2006); MellemaeLal] (2006b) and Il iev et alJ 
(12007 al) . While this volume is too small to allow us to consider the 
rarest , most luminous sources like the SDSS QSO's jpan & et alJ 
I2OOI l2006h . we have sufficient resolution to locate the majority 
of sources responsible for reionization and take explicit account of 

their radiation, and to derive good quality absor ption spect ra. 

Of the range of simulations presented in llliev et alJ (l2007ah 
we here consider one specific run, labelled f250C. Our simulations 
were performed using a combination of two very efficient compu- 
tational tools, a c osmological particle-mesh code called PMFAST 
dMerz et al.l2005h for following the structure formation, whose out- 
puts are then post-processed using our radi ative transfer and non - 
equilibrium chemistry code called C^-Ray j Mellema et alj|2006ab . 
The parameter = 250 characterizes the emissivity of the ion- 
izing sources - how many ionizing photons per gas atom in the 
(resolved) halos are produced and manage to escape from the host 
halo within ^ 20 Myr, which is the time between two consecu- 
tive density slices, equal to two radiative transfer timesteps, while 
'C indicates that this run models the gas clumping at small (sub- 
radiative transfer grid) scales based on a fit given by 

^ / \ o/?oni^ -0. 1822Z+0. 003505 ^-.x 

Csuh-gTid{z) = 26.2917e ^ . (1) 

based on the used WMAP3 cosmology (a good fit for 6 < z < 30). 
This fit to the small-scale cl umping fac tor is a more precise version 
of the one we presented in llliev et aP ( 2005). We derived it based 
on a PMFAST simulation with a computational mesh of 3248^ 
and with particle number of 1624^, with computational volume of 
(3.5 Mpc)^. These parameters correspond to particle mass of 
10^ Mq, minimum resolved halo mass of lO^M©, and a spatial 
resolution of 1 kpc comoving. 

Our (100 Mpc)^ volume simulation resolves all halos 
with mass above 2.2 x 10^ M©, higher than the mass above which 
atomic-line cooling of hydrogen becomes effective, which is ^ 
10^ Mq . As a consequence, our treatment does not include the con- 
tribution of low-mass sources to reionization. Higher-resolution, 
smaller-box simulations which do i nclude all ionizing sources 
above the atomic-line cooling limit (llliev et al.l 12007 ah indicate 
that the effects of low-mass sources are primarily confined to the 
earliest stages of reionization, when such sources are dominant. 
Throughout most of the reionization process, and especially during 
its late stages, the low-mass sources, which are strongly clustered 
around the high density peaks, are heavily suppressed due to Jeans- 
mass filtering in the ionized regions and thus have limited effect on 
the reionization progress and large-scale geometry. Since the hy-a 
observations largely probe the later stages of reionization, where 
the n eutral gas fraction is ~ 30% or less jMalhotra & Rhoadl 
I2OO4"), we do not expect that our conclusions will be strongly af- 
fected by the abse nce of low-mass sources. Larger simul ations, cur- 
rently in progress ( Iliev et al.l2008bl : IShapiro et al.l2008 ). which re- 
solve all atomically-cooling halos in ~ 100 Mpc boxes will settle 
these uncertainties. 

Throughout this work we assume a flat (Qfc = 
0) ACDM cosmology (^^m, ^^a, ^^b, /i, erg, n) = 
(0.24, .76, 0.042, 0.73 0.74 , 0.95) based on WMAP 3-year 
results ( Spergel & et"ani2007h . hereafter WMAP3. Here Qm, ^a, 
and Qb are the total matter, vacuum, and baryonic densities in 
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Figure 1. Mass accretion history of the three most massive halos found in 
our computational volume at redshift z = 6. The mass growth is roughly 
exponential in redshift, and is well-approximated by M(z)[Mq] = 
exp(A — az), where A = 31.2, a = 0.52 for most massive halo (solid, 
red), A = 32.5, o; = 0.81 for the next most massive halo (short-dashed, 
blue) and A = 32.0, o; = 0.75 for the third most massive halo (long- 
dashed, green). Fits are shown by the thin straight lines (with corresponding 
line types and colors). 



the linear perturbation theory, and n is the index of the primordial 
power spectrum of density fluctuations. 



3 LUMINOUS HIGH-REDSHIFT SOURCES AND THEIR 
ENVIRONMENT: PROPERTIES, EVOLUTION AND 
REIONIZATION HISTORY 

The luminous sources at high redshift, the ones typically seen in 
current surveys, are hosted by rare, massive halos which form at the 
location of the highest peaks of the density field. The statistics of 
Gaussian fields predicts that such high density peaks are rare and 
highly clustered in space, more strongly so at high redshifts. As 
a consequence, each high-redshift, massive galaxy should be sur- 
rounded by numerous smaller ionizing sources. The self-consistent 
reionization history simulations of such regions require following 
a sufficiently large volume, in order to obtain the correct statis- 
tics and biasing of the rare peaks, while at the same time resolving 
all the low-mass halos which are the main drivers of the reioniza- 
tion process. Our current radiative transfer simulations are able to 
achieve this. We also note that correctly modelling the nonlinear 
bias of the rare peaks in semi-analytical modes is a difficult and still 
unsolved problem. As a consequence, the halo clustering in semi- 
analytical models is typically underestimated, and in some cases 
even ignored. This often yields incorrect results, e.g. in estimates 
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3.1 Mass accretion history of massive halos at high redshift 

In Figure [T] we show the mass-growth history of the three most 
massive (at redshift z = 6) halos found in our computational vol- 
ume, with final masses of 1.5 x lO^^M©, 9 x lO^^M© and 8.2 x 
10^^ M0, respectively. All correspond to very rare, ^ 4.5 — 5 — a 
peaks of the density field. The first progenitors of these halos re- 
solved in our simulation (Mhaio > 2 x lO^M©) form very early, 
at z ^ 16. Thereafter, the halo masses grow roughly exponen- 
tially with redshift, M oc exp(— az). This behaviour is in good 
agreement with previous results on the growth of dark m atter ha- 
los in hierarchical ACDM models dWechsler et al.ll2002h . as well 
as wi th similar results o btained in more generic, idealized situa- 
tions ("Shapiro et al .l2004h . The slopes a for our halos exhibit much 
smaller scatter (0.52 < a < 0.81; see Figure [T] caption for the 
complete fitting formulae) than the ones found by Iwechsler et al.l 
(120021) (0.4 < a < 1.6). The exponential fits are excellent during 
periods when no major mergers (i.e. mergers with halos of simi- 
lar mass) occur. This is the case e.g. at late times (z < 11) by 
which time these halos have grown enough to dominate their sur- 
roundings and no similar-mass halos remain in their vicinity. The 
exception here is the second most massive halo (short-dashed line). 
Its last major merger occurs ai z ^ S. At early times (z > 12) the 
mass evolution of all three halos does not follow the overall expo- 
nential trend. Instead, the most massive halo is growing faster than 
its long-term trend, through series of rapid major mergers, while 
the other two halos initially grow slower than their long-term trend. 



3.2 Reionization history 

In Figure [2 (top and third row) we illustrate the main stages of the 
reionization history of a high density peak and its intergalactic envi- 
ronment. The particular peak we show here is the one surrounding 
the largest-mass halo found in our computational volume at red- 
shift z — 6. For clarity and convenience for the reader we have 
shifted the peak to the box center using the periodicity of our com- 
putational volume. The first resolved source in this region forms 
ai z ~ 16. By redshift z = 12.9 (top left; mass- weighted ionized 
fraction is Xm = 0.001 at this time) the central halo has already un- 
dergone several major mergers with nearby halos and its total mass 
has grown to 1.5 x 10^° M©. The combined effect of the central 
halo and the nearby members of the same source cluster is to cre- 
ate a substantial ionized region (of size R ^ 1 h~^Mpc, defined 
as the radius at which the integrated continuum optical depth from 
the source at the Lyman limit reaches unity). At this time there 
are only a few ionized regions in our computational volume (and 
just two intersect the plane shown). By z = 10.1 (top middle; 
Xm — 0.10) many more H II regions appear and both the sources 
and the ionizing regions are strongly clustered in space. The cen- 
tral region still remains the largest one. By redshift 2; = 9.0 (top 
right; Xm — 0.28) many more haloes have formed, most of them in 
large clustered groups. The H II region surrounding the central peak 
remains among the largest {R ^ Q h~^Mpc), but several other ion- 
ized bubbles reach comparable sizes. 

At redshift z — 7.9 (bottom left; Xm — 0.66) some quite 
sizable regions, more than ten Mpc across, have percolated, reach- 
ing local overlap. The central bubble has grown to a size of R ^ 
8h~^Mpc. The reionization geometry becomes quite complex, 
with most ionized bubbles becoming interconnected, while leav- 
ing large neutral patches in-between. By 2; = 7.0 (bottom mid- 
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Figure 3. Histogram of the bubble size distributions along LOS. Shown 
are the distances from the source (the most massive halo in our volume 
at 2; = 6) at which continuum optical depth (at the hydrogen ionization 
threshold) surpasses r = 1 (solid lines), or r = 4.6 (dotted lines), at 
several illustrative redshifts, as labelled. 



already merged into one topologically-connected region, although 
substantial neutral patches still remain interspersed throughout our 
volume. The volume remains on average quite optically-thick to 
ionizing continuum radiation. Finally, at 2; = 6.0 (bottom right; 
Xm = 0.9999; R = 18.7h"^Mpc) our volume is well beyond 
overlap (which we define by Xm > 99%). Only by that time the 
volume becomes on average optically-thin to ionizing radiation. 

The corresponding panels in the second and fourth row of Fig- 
ure |2] show images of the (non-equilibrium) ionization rate distri- 
bution at the same redshifts. The distribution is highly inhomoge- 
neous, following the patchiness and peaking strongly in the vicin- 
ity of large source clusters. The volume- averaged photoionization 
rates at redshifts z = (12.9; 10.1; 9.0; 7.9; 7.0; 6.0) in units of 
10"^^ s-^ are r_i2 = (7.0 x 10"^; 1.3 x 10"^ 2.3 x 10"^ 7.6 x 
10~^;1.3;4.9), growing strongly with time as larger fraction of 
the volume becomes ionized and ever more sources form. Except 
for the earliest times, the peak photoionization rate values on the 
grid (corresponding to the brightest points of the images) remain 
fairly constant with time, at the same redshifts they are r_i2 = 
(1.0 X 10^; 2.0 X 10^ 2.0 x 10^; 2.5 x 10^; 2.8 x 10^ 5.1 x 10^). 
The photoionization rate distributions and evolution are discussed 
further in ^ 1441 

To further characterize the shape and boundary sharpness of 
the central H II region in Figure[3]we show a histogram of the radial 
distance R from the central source at which the cumulative contin- 
uum optical depth at the ionizing threshold of hydrogen along each 
LOS reaches unity (black) and 4.6 (red; this optical depth value cor- 
responds to 1% transmission). A spherically- symmetric H II region 
would correspond to a single value for the radius for a given optical 
depth, regardless of the directionality, while any scatter around the 
peak would be a measure of the non- sphericity of the ionized re- 
gion. Furthermore, a comparison between the distributions for the 
two optical depths measures the sharpness of the H II region bound- 
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Figure 2. The reionization history of a high density peak. The images are centered on the most massive (at z = 6) halo in our computational volume and are 
of size 100/z~^Mpc to the side. The snapshots are at: (top row) z = 12.9 (left; global mass-weighted ionized fraction Xm = 0.001), z = 10.1 (middle; 
Xm = 0.10), z = 9.0 (right; Xm = 0.28), and (third row) z = 7.9 (left; Xm = 0.66), z = 7.0 (middle; Xm = 0.94), and z = 6.0 (right; Xm = 0.9999). 
The underlying cosmological density field (dark,green) is superimposed with the ionized fraction (light, orange) and the cells containing ionizing sources (dark, 
blue dots), the slices thickness of 0.5 Mpc (1 cell) in the density and ionized fraction fields and 10 Mpc in terms of sources. The corresponding 
images of the (non-equilibrium) photoionization rates (0.5 Mpc thickness) are shown in the second and bottom rows. 
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Figure 4. Projected distributions at z = 7 of halos within 1 h -"^Mpc comoving (left) and within 5 h -"^Mpc comoving from the most massive object in the 
box. Circle areas are proportional to the mass of the corresponding halo. 



passed simultaneously (since the neutral gas is extremely optically- 
thick) and in such situation the histograms would coincide. On the 
other hand, if the boundary of the ionized region is not clearly de- 
fined due to local percolation, gas is highly ionized and the optical 
depth along the LOS increases only slowly, thus the values of 1 and 
4.6 are reached at very different distances from the source. 

At redshift z — 9 (Xm — 0.28)the H 11 region surrounding the 
central source is largely spherical and its boundary is well-defined, 
albeit with some modest spread around the peak, indicating some 
departures from sphericity, in agreement with what was seen in 
Fig.ia At redshift z = 8 (xm = 0.62) the H 11 region remains 
fairly spherical, slightly more so for r = 1. The r = 4.6 histogram 
has a long tail at large values of R, up to ^ 25h~^ Mpc, meaning 
that a small percentage of the LOS reach 99% opacity only at such 
fairly large distances. Both distributions are somewhat wider than 
before and start to depart from each other, reflecting the increasing 
"fuzziness" of the bubble boundary as it merges with other nearby 
bubbles due to local percolation. 

At z = 7 (Xm = 0.94) the T = 1 distribution still re- 
tains a well-defined peak, albeit one accompanied by a long high- 
R tail reaching out to tens of Mpc. However, the r = 4.6 his- 
togram changes its character completely, becoming very broad, 
with only a low peak close to the r = 1 peak and several sec- 
ondary peaks. A significant number 14%) of the LOS do not 
reach r — 4.6 within 50h~^Mpc of the central source (these are 
collected in the last, R — 50h~^Mpc bin of the histogram). Fi- 
nally, at 2; = 6 {xm = 0.9999), well beyond the overlap epoch 
{z ^ 6.6, Xm = 0.99) there is no indication of a defined ionized 
bubble. Both distributions are very broad, with no clear peak. At 
this time the optical depth is determined by the tiny remaining neu- 
tral fraction, which in its turn is dictated by the density variations 
of the Cosmic Web. Either value of r is reached at a wide range 
of distances in different directions, from ^ 10h~^Mpc to over 
50h~^Mpc. The IGM becomes largely optically-thin to ionizing 
radiation and thus for most LOS neither t — 1 nor r = 4.6 are 
reached within 50 h~^Mpc from the source. 




R [h-iMpc] 

Figure 5. Histograms in radial bins of the total cumulative emissivity from 
all halos within distance R, A^ph,tot(^)' measured from the central ob- 
ject, here the most massive object in the computational volume, in units 
of the central object's emissivity, A^ph, central- Shown are (bottom to top 
on the right) z = 12.9, 9, 8, 7 and 6.25 (xm = 0.001, 0.28, 0.62, 0.94 
and 0.998). The last radial bin contains the total emissivity of the computa- 
tional box at that time. The vertical lines on top (corresponding line types 
and colors) indicate the average size of the H II region at the corresponding 
redshift. 



hard sources, e.g. QSOs, contribute significantly to reionization the 
ionized region boundaries will inevitably be thicker and the transi- 
tion from ionized to neutral smoother, due to the longer mean free 
paths of the hard photons. However, most current observational in- 
dications are that stellar sources dominate reionization, in which 
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3.3 Halo clustering in the vicinity of a luminous source 



Taking a closer look at the halo clustering in the immediate vicinity 
of a high density peak in Fig.|4]we show the projected spatial distri- 
bution of halos around the most massive halo at z = 7. There are 31 
resolved halos within 1 /i~^Mpc from the most massive halo and 
360 resolved halos within 5 h~^Mpc (both including the largest 
halo itself). The area of each circle is proportional to the mass (and 
thus, in accordance to our source model, to the luminosity) of the 
corresponding halo. Halos are distributed very anisotropically, con- 
centrating preferentially along the density filaments and sheets of 
the local Cosmic Web. 

The mass of the most massive halo is 9 x 10^^ at that time, 
well above any other halo in its vicinity. Nonetheless, the low-mass, 
but numerous halos surrounding the peak contribute significantly to 
the total ionizing emissivity. In order to quantify this point further, 
in Figure |5] we show the cumulative emissivity vs. radial distance 
from the central halo for several redshifts spanning the whole range 
of interest here. At all redshifts only within ^ 2 h~^Mpc from it- 
self does the most massive halo dominate the total photoionizing 
emission (i.e. it contributes more than 50% of the cumulative to- 
tal emission coming from that region). The exception is the high- 
est redshift (z — 12.9), in which case the surrounding cluster of 
sources, rather than the central source, dominate the total flux even 
within 1 h~^Mpc of the peak. The reason for this is the presence 
nearby of another halo of almost the same mass as the central one, 
which shortly thereafter merges with it. The cumulative emission 
is dominated by the small halo contribution beyond that distance. 
Within 10h~^Mpc the central source contributes only 10-30% of 
the total emission. As more and more halos form and the lowest 
mass ones become relatively common the fractional contribution 
of the most luminous source to the total emissivity gradually de- 
creases in time for all radii larger than a few comoving Mpc. At 
30h~^Mpc the central source is dominated by the rest by 1.5- 
2 orders of magnitude. As fraction of the total emissivity of all 
sources in our computational volume (shown in the last bin of the 
histogram) the most luminous one contributes only ~ 2% of the 
total at z = 12.9, decreasing to 0.1% at 2; = 6.25. 

In Figure [5] we also indicate the current H II region mean ra- 
dius (vertical lines of corresponding line types and colors). Within 
its own bubble the central source contributes ~ 50% of the emis- 
sion at z = 12.9, decreasing to ~ 10% at 2; = 7 and ^ 1% at 
z = 6.25. 



Recently, IWvithe & Loebl (l2005h presented a semi-analytical 
model of the source clustering at high redshift. They found that a 
central galaxy at z = 7 and with a velocity dispersion similar to our 
most massive source, av ~ 200 kms~^) contributes ^ 40 — 70% 
of the H II region radius, or ^ 10 — 35% of its volume (see their 
Fig. 2, right panel), a factor of a few larger compared to our results 
where we find that the central galaxy contribution to the total flux 
is ~ 1 — 10%. This discrepancy is most probably due to underes- 
timate of source clustering in their bias model, and possibly also to 
the differences in the assumed source efficiencies. Furthermore, we 
have only considered a single massive source, which of course does 
not account for random statistical variations from source to source. 
Nevertheless, this result underlines the importance of considering 
the luminous sources individually (as opposed to considering an 
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Figure 6. Spherically-averaged profiles of (bottom to top on the left/right): 
the comoving number density of hydrogen (in cm~^; short-dashed, 
blue), the neutral hydrogen fraction (solid, red), radial velocity, Vr (in 
100 km s"-*^; long-dashed, green) and continuum optical depth at hv = 
13.6 eV integrated from the source (R = 0) outward, all at redshift z = 6.6 
{xm = 0.99) vs. R, the comoving radial distance from the most massive 
galaxy. 
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Figure 7. Mean radial velocity at 2; = 6 of the IGM with respect to the 
source (red solid line, negative means towards the halo, positive - away) 
and its rms variation (error bars), both plotted vs. (comoving) radius from 
the source. 



3.4 IGM Environment of luminous sources at high redshift 

The Ly-a absorption is strongly influenced by the local IGM en- 
vironment - density, velocity and ionization structure - around the 
the source. As an illustrative example, let us consider a particular 
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simulation is already highly-ionized, with a mean mass-weighted 
ionized fraction of Xm = 0.939, and a volume ionized fraction of 
Xv = 0.925. In Figure[6]we show the spherically- averaged profiles 
around the central source of the gas number density, n, neutral frac- 
tion, xhi, integrated continuum optical depth (r, at hu = 13.6 eV) 
and radial velocity, Vr. The radiative transfer cell which contains 
the central source is highly-overdense, with 6 = n/n — 1 = 256, 
which indicates that this cell roughly coincides with the source 
halo itself. The radial density profile declines steeply away from 
the source. The overdensity is ^ = 3.3 slI R ^ lh~^Mpc, de- 
creasing to^ = lati?^2.5 h~^Mpc, and approaching the mean 
density at distances beyond R ^ 10 h~^Mpc. 

The radial velocity profile shows an extended (R ^ 
20h~^Mpc) infall region, with the mean radial velocity peaking 
at ^ 150 kms~^ before dropping to zero inside the source halo it- 
self. The proximity region of the central source (R < 15 h~^Mpc) 
is highly-ionized, with neutral fraction xhi < 10~^. The rest of 
the volume still has an appreciable neutral fraction 0.1 — 1%), 
however, and is thus on average still optically-thick, with the mean 
optical depth reaching r = 63 at = 50 h~^Mpc. 

However, the spherically-averaged quantities provide only a 
limited information about the state of the IGM surrounding each 
source. All quantities are distributed highly-anisotropically, and 
thus affect the Ly-a emission differently along each LOS. In par- 
ticularly, the effect of the relative velocities of the IGM and source 
is relatively poorly studied at present. An optically-thick medium 
at rest with respect to a Ly-a source would absorb the blue wing of 
the line and transmit its red wing , at longer wavelengt hs than the 
line centre at Ao = 1215 A (e.g. Miralda-Escuddl 1 998b . A relative 
motion of the IGM gas and the source along the LOS would result 
in either more or less transmission, depending on the motion di- 
rection. E.g. gas infall towards the source along the LOS from the 
observer would redshift it, resulting in some absorption of the red 
wing of the line. We note that in order to evaluate these velocity 
effects with any precision much higher resolution simulations (and 
ones including outflows) will be required. Our radiative transfer 
grid has cell size of ^ 0.5/ h Mpc comoving, or ^ 70//i kpc phys- 
ical at z = 6, which roughly corresponds to the size of the largest 
halos found in our box. The velocity and density fields used are at 
twice higher resolution 0.25/ h Mpc comoving). Therefore, our 
results here should be considered as a guidance, illustrating that the 
velocity effects are quite important and should not be ignored. 

In Figure [7] we show the average IGM gas velocity relative 
to the source (line) vs. distance from it and the variance of that 
average velocity (error bars). As noted above, in the vicinity of the 
source the IGM on average infalls towards the halo. However there 
are large variations, of order hundreds of km around this mean. 
E.g. a velocity offset of 200 km at z = 6 corresponds to A A ^ 
6 A, of the same order as the typical observed line widths, and thus 
a relative motion of the IGM and source of this order could have a 
very significant effect on the observed line. In the next section we 
would quantify the effect of peculiar velocities on the Ly-a line. 



4 OBSERVABILITY OF HIGH-Z LY-a SOURCES 

4.1 Absorption spectra of luminous sources 

Much of the information about high-redshift Ly-a sources and IGM 
is based on absorption spectra. We thus start our discussion of the 



our computational volume. Figures [8l[TT] we show sample spectra 
along three random lines-of- sight at a few selected redshifts span- 
ning the complete range of interest here. On the left panels we show 
the distributions of Ly-a Gunn-Peterson optical depth, tgp, neu- 
tral fraction, xhi = 1 — x (multiplied by 10^ for clarity), and 
gas density in units of the mean, A = n/n. On the right panels 
we show the corresponding Gunn-Peterson transmission spectra for 
flux level of unity, exp(— tgp). The horizontal lines on the left pan- 
els indicate the optical depth of r = 4.6, which corresponds to 1% 
transmission. For reference, this value is roughly equal to the opti- 
cal depth of a hydrogen gas with neutral fraction of xhi = 10 ~^ 
at the mean density at redshift z = 6.6. All quantities shown are 
in redshift/wavelength space and in the observer (z = 0) frame. 
For reference, on the top axis of each figure we show the approx- 
imate corresponding distances in real space. Finally, in Figure [TJl 
we show the mean, averaged over all random LOS, transmission 
spectra at the same redshifts, for the most massive source (left) and 
mean over all sources (right). 

The Lyman-a absorption as a func tion of wavelength i s com- 
puted using the standard procedure (e.g., iTheuns et al .l (119981)). The 
optical depth and transmission results include the redshift- space 
distortions due to the local peculiar velocities, relative to the pecu- 
liar velocity of the source (i.e., after applying peculiar velocity dis- 
tortions, the whole spectrum has been shifted slightly so the source 
is returned to it's real space position). The temperature of the gas is 
assumed to be 10^ K when computing thermal broadening, consis- 
tent with the assumption adopted for the simulations. 

The nominal resolution of our spectra is R ^ 6000 — 12000 
at z = 6 (higher at higher redshifts), based on our grid reso- 
lution of 203^ (radiative transfer) and (406^ density and veloc- 
ity fields). This resolution roughly corresponding to the one for 
medium-resolution observed spectra. In reality the situation is more 
complicated. Due to the non-linear transformations between our 
raw simulation data and the final spectra the limited simulation res- 
olution can affect the results even if it were better than the obser- 
vational resolution. A separate issue pointing in the same direction 
is the fact that the real data has effectively infinite resolution in 
the transverse direction (i.e., the width of the light beam), and thus 
in order for us to be accurate we have to ensure we are resolving 
essentially all the transverse structure, which is not the case for 
the current simulations. As a result, our spectra should not be con- 
sidered completely realistic predictions, but rather as a guidance 
showing some important features to be expected from real spectra, 
as discussed below. Future higher-resolution, more detailed simula- 
tions will be better suited to make realistic predictions of the actual 
detailed spectral properties. 

At early times (z = 12.9; Figure [8j top) the H II region 
surrounding the most massive source is still quite small (see also 
Figs.[2l[3]and[5]) and the Ly-a emission of the source is completely 
suppressed by the damping wing of the Ly-a line profile, rendering 
it unobservable. The ionized region grows quickly after that and by 
z = 9 reaches size of 6 h~^Mpc (Figure[8j bottom). Regardless, 
the Ly-a optical depth even within the source proximity region re- 
mains quite high, at few up to ~ 10, which allows through only a 
very weak transmission. The damping wing slightly weakens com- 
pared to the higher redshifts, but is still quite substantial, and still 
depresses most of the red wing of the emission line. 

Some of the continuum immediately behind the source in red- 
shift space (within a few A) is absorbed due to gas infall towards 
the density peak. Because of that additional velocity towards the 
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Figure 8. Sample LOS at redshifts z = 12.9 (top; Xm = 0.001) and 9.0 (bottom; Xm = 0.28) vs. A/comoving distance from the most massive galaxy. 
Shown are (left panels) the optical depth (solid), neutral fraction xui = 1 — x (xlO^; dotted) and density in units of the mean (dashed), and (right panels) 
the corresponding transmission. The vertical lines show the position of the central source (in redshift space, i.e. accounting for its peculiar velocity along the 
LOS). The horizontal lines on the left indicate the optical depth equivalent to 1% transmission. On the right, the shaded region is the transmission in the case 
where the unabsorbed spectrum is flat (the horizontal dotted line). 



the luminous source's position is completely dark due to the high 
density in the middle of the density peak, regardless of the very low 
neutral fraction there. These features persist throughout the evolu- 
tion and are very characteristic for all luminous sources, since these 
always associated with high density peaks and surrounded by infall. 

At redshift z = 8.1 (Figure |9] top) a weak damping wing is 
still present and essentially no transmission occurs on the blue side 
of the line. The sources at this time may be potentially visible with 
very deep observations. Only by redshift z = 7 (Figures [12]) the 
ionized region is sufficiently large for the damping wing to effec- 



photons on the blue side of the Ly-a line along most LOS. On 
average a weak transmission at a few percent level starts coming 
through in the proximity region of the source (Figur^TJ}. As more 
and more sources form and the ionizing flux rises the first trans- 
mission gaps start to appear at 2; < 7, both in the mean IGM away 
from the peak and in the source proximity (Figure [TO]). At redshift 
z = 6.6, our nominal overlap time (Figures [TO] and (TJ]), the prox- 
imity region extends for 30 A and has become fairly optically- 
thin, allowing up to 30-40% transmission. Most of the volume still 
remains optically-thick, but some substantial transmission regions 
appear in the IGM away from the peak. We also note that there are 
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ing for much more transmission than others. The size and prop- 
erties of the proximity region also vary due to its asymmetry and 
the anisotropy of nearby structures. Finally, during the post-overlap 
epoch (Figures [Tol bottom, [TT] and [12]) the IGM slowly becomes 
more optically-thin to Ly-a and gradually approaches the state of 
the Ly-a forest. There are no more clearly-defined H II bubbles. 
Only a few isolated low-density regions remain neutral. This is due 
to the inside-out character of the reionization process, whereby the 
high-density regions are preferentially ionized first, while the voids, 
where structure formation is delayed are ionized last. 

Comparing the radially- averaged mean transmission around a 
luminous source (high density peak) and the average for all sources 
(i.e. the mean behaviour around a typical source; Figure [12]) we 
see both some similarities and several notable differences. The 



(z ^ 8 — 9) and finally disappearing at later times (z < 7). Nai- 
ively, one might expect that compared to an average source the lu- 
minous ones would suffer from weaker damping since such sources 
typically reside in the middle of large H II regions, far from their 
boundaries. In fact, this expectation proves only partially correct. 
The damping is somewhat more pronounced around an average 
source than around a luminous one, but the differences we ob- 
serve are rather modest, reflecting the fact that weaker sources are 
strongly clustered around the same density peaks that host the lu- 
minous ones, thus largely share the damping (or its lack) with the 
central source. For the same reason the damping becomes irrele- 
vant at about the same time in both cases, which would not have 
been the case if the weak sources were residing in smaller, isolated 
bubbles. 
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of the source. The same behaviour is not seen for a typical source. 
Such smaller halos tend to move more in tandem with their sur- 
rounding IGM, often towards the nearest high density peak. Some 
local infall should exist also for these halos, but this is at very small 
scales, unresolved here. However, these scales are small compared 
to the typical emission line width (see next section) and thus we do 
not expect that such local infall has significant effect on the emis- 
sion line. 

One final important difference between a luminous and an av- 
erage source is that the latter does not typically have a proxim- 
ity transmission region on the blue side of the line. The spectra 
of the luminous sources, on the other hand exhibit extended high- 
transmission (10-60% transmission) regions within 5 Mpch~^ (^ 
30 A). This behaviour is again due to the high source clustering 



gas density there. The line center coinciding with a high density 
peak remains optically-thick, however, unlike the line center of a 
typical source. Away from the proximity region the absorption is 
largely saturated, but there are a number of transmission gaps with 
up to a few per cent transmission. Future work would quantify the 
statistics of these features and its evolution. 



4.2 Emission line shape and its evolution 

In order to study the effect of IGM absorption on the line profile 
shape we model the intrinsic Ly-a line as a Gaussian with an rms 
width of 160kms~^ and peak amplitude normalized to unity. In 
Figures (TS] [14] and [15] we show sample results for several LOS 
through the most luminous source at redshifts z = 8, 6.6 and 6, re- 
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Figure 11. Same as Fig.[8] but at redshift z = 6.0 (xm = 0.9999). 
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Figure 12. Mean radially-averaged transmission around the most massive source (left) and average over all sources (right) at several representative redshifts, 
as labelled. The corresponding ionized fractions by mass are Xm = 0.001, 0.105, 0.279, 0.618, 0.939, 0.992, 0.9986 and 0.9999. 



sumed intrinsic (black) transmitted (red) emission line, while on 
the right panels we show the corresponding distributions of Ly-a 
Gunn-Peterson optical depth, tgp, neutral fraction, xhi = 1 — x 
(multiplied by 10^ for clarity), and gas density in units of the mean, 
A — n/n and again (for reference) the intrinsic emission line. 

At redshift z = 8.1 (Figure [13]) the emission line shape is 
fairly regular and does not vary strongly for different LOS. The blue 
wing is generally highly absorbed and no appreciable flux comes 
through and much of the red wing is absorbed as well. The reasons 
for this behaviour become clear from the right panels. The neutral 



high density of the peak and its immediate vicinity results in optical 
depths of TGP > 10 everywhere and r > 100 at the peak itself. The 
gas infall towards the peak leads to a significant absorption of the 
red wing. The only transmitted flux comes from the far red side of 
the line (slightly depressed by the weak remaining damping wing). 

At overlap (1% global neutral fraction by mass; z = 6.6) the 
line shape becomes much more irregular and varies significantly 
between the different LOS (Figure[T4]). Significantly larger fraction 
of the flux is transmitted, both on the red and on the blue side of 
the line. The neutral fraction in the vicinity of the source is still 
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Figure 13. Emission lines: (left panels) intrinsic line, assumed a Gaussian with rms width of 160 km (black, top), and transmitted one (red, bottom) for 
three sample LOS at redshift 2; = 8.1 (xm = 0.62), and (right panels) the corresponding optical depth (solid, black), neutral fraction xui = 1 — x (xlO^; 
dotted, red) and density in units of the mean (dashed, green), all in redshift space (i.e. accounting for the relative velocities). The intrinsic emission line is also 
shown for reference. 
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Figure 14. Same as Fig.[T3] but at redshift z = 6.6 {xm = 0.99). 



distribution. The effect of the gas infall towards the peak is still 
present and some of the red wing is absorbed, but much less so than 
at higher redshifts since the infalling gas is more highly-ionized. 
The line center remains absorbed and all damping wing effects have 
disappeared. 

After overlap (Figure [151) the neutral fraction gradually de- 
clines, decreasing the optical depth and allowing ever more flux to 
be transmitted. The neutral fraction is still fairly uniform and thus 



the different LOS. There is transmission in both red and blue wing 
of the line. 

In Figure[T6lwe show the evolution of the mean (i.e. averaged 
over all LOS) observed emission line shape for the most luminous 
source in our volume (left) and average over all sources (right). 
In both cases the line starts completely damped {z — 12.9). The 
later evolution of the mean line shape differs significantly, how- 
ever. By redshifts 2; = 9 — 10 a significant fraction of the red wing 
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Figure 16. Evolution of the mean emission lines for most massive source (left panels) and average over all sources (right panels). Shown are the intrinsic 
emission line (dotted, red; assumed a Gaussian with rms width of 160kms~^, normalized to one at the peak), the transmitted line (solid, red), and mean 
absorption (solid, black) for three sample LOS several redshifts, as labelled. 



average source. At later times (z 7 — 8) this situation is reversed 
- practically all of the red wing of the line is transmitted for an 
average source, but much of the flux is still absorbed for the lu- 
minous source due to the high density peak in the middle and its 
surrounding infall. As a word of caution we should note that some 
of this effect is in fact not physical but numerical, since our simu- 
lations do not resolve well the detailed structure around the smaller 
halos. However, as we also mentioned above, these resolution ef- 
fects should be modest considering that the emission line is fairly 
wide and thus reasonably well-resolved and any corrections due to 
smaller-scale structures will not affect much of the line. 



ences between the luminous and average sources. The strong clus- 
tering of sources around the density peaks result in very high fluxes 
and thus a more pronounced highly-ionized proximity region blue- 
ward of the line center. Thus, significantly more of the blue wing 
of the luminous source line is transmitted, up to 10% on average at 
2; = 6, vs. only ^ 1 — 2% for an average source. 

An interesting consequence of the very high absorption ob- 
served at the line center for massive sources and the redshift- space 
distortions due to the gas infall (the latte r similar to the one studie d 
theoretically in a more idealized setup bv lBarkana & Loebl (l2004ah ) 
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Figure 17. Overall mean transmission fraction (not necessarily near a 
source) in Lyman- o; (solid, black), Lyman-/? (dashed, red), and Lyman-7 
(dotted, g reen). The points with horizontal error-bars show the measure- 
ments of (iFan & et'^l2006h . The (black, red), (lower, upper), points with 
vertical error bars show Ly-(a,(3), while the green dot shows Ly-7 (see text). 



Figure 18. PDF of the photoionization rate at at z = 10.1 (blue; Xm = 
0.105), z = 7.0 (red; Xm = 0.94), and z = 6.0 (green; Xm = 0.9999). 
We show the actual, non-equilibrium rates (solid) and the corresponding 
equilibrium rates (dotted, same color at each redshift). All PDF's are nor- 
malized to have an area of unity below the curve. 



in principle it might be possible to use the line profiles of bright Ly- 
a sources to study the infall surrounding their host halos. In prac- 
tice this might be difficult due to a number of complications. The 
line structure is quite different along different LOS, partly (as we 
pointed above in § 2.4) due to the very anisotropic velocity structure 
surrounding the source, as well as its own peculiar motion. Further- 
more, our analysis only takes into account the effects of the IGM on 
the line shape, while realistic ones will be affected also by the host 
galaxy's internal structure, outflows, etc. Modelling all those ef- 
fects correctly will require high-resolution radiative-hydrodynamic 
simulations, which is well beyond the scope of this work. 



4.3 Evolution of the mean transmissivity 

Figure [T7] shows the mean transmission fraction as a function of 
redshift, for the Lyman-a, /3, and 7 transitions. The latter two are, 
respectively, 6.2 and 17.9 times weaker than Lyman-a, so trans- 
mission can be see n in cases where Lyman-alpha would be opaque 
(iFan & et al.l l2006^. At the low redshift end of our simulation data 
these q uantities have been observed in the hig hest redshift quasar 
spectra (iFan & et al.|[2006h . Fig. [17] shows the (iFan & et alJl2006h 
(points with error bars). At 2; > 6, the Ly-a and Ly-/3 transmis- 
sion measurements are only upper limits. It is unclear exactly what 
error bar should be assigned to the Ly-7 point. Fan et al. (2006) 
presented two upper limits and a detection in Ly-7, which, taken 
together, support a detected mean transmission level (plotted in the 
figure) well below our prediction, with relatively small (^ 40%) 
errors; however, with only these three points we can not be sure 
that the sample variance isn't substantially larger. It appears that 
the reionization model in our simulations reproduces roughly the 
correct tail-end of reionization, but the data favor somewhat less 



4.4 Photoionization Rates 

In Figure [18] we show the normalized probability density distri- 
butions (PDFs) of the nonequilibrium photoionization rates for all 
cells in our computational volume at three representative redshifts - 
z 10.1 (xm 0.105), 7.0 (xm 0.94) and 6.0 (xm 0.9999) 
(early times, late times and well after overlap) in linear (top) and 
log (bottom) scales. For comparison we also plot the photoioniza- 
tion rates if the corresponding cells were in ionization equilibrium. 

The right peak of each PDF distribution reflects the most com- 
mon photoionization rate values in the ionized regions. The peak 
position remains at r_i2 ^ 1 throughout the evolution, with only 
slight shifts. The distributions are quasi-Gaussian, but with long 
non-Gaussian tails at both high and (especially) low values of F. 
As could have been expected, the fraction of high-F cells, all of 
which either contain sources or are in the immediate vicinity of a 
source, grows strongly with time, as ever . The highest photoioniza- 
tion rate values we find reach r-12 10^ — 10^. This peak value 
rises over time, as a consequence of the growth of galaxies and the 
large number of sources forming in and around the density peaks. 
In the ionized regions the equilibration time is short and photoion- 
ization equilibrium is generally a good approximation. The cells 
with lower values of F (below r_i2 ~ 0.01 — 0.1 correspond to 
the ionization fronts or neutral regions. There are many more cells 
in I-fronts at z = 10.1 than at later times, when most of the IGM 
is already ionized. The equilibrium and actual rates differ widely in 
those regions, indicating that the assumption of ionization equilib- 
rium would be a very poor approximation there. 

The photoionization rate-density correlations at z = 9 and 
2; = 6 are shown in Figure [19] as contour plots. At high densi- 
ties there is a clear, and fairly tight, correlation between the density 
and the photoionization rate. This reflects the fact that many more 
sources form in overdense regions, resulting in higher photoioniza- 
tion rates. At densities just above the mean the correlation becomes 
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Figure 20. (bottom panels) Ly-a luminosity function of high-redshift sources without (black) and with absorption included (red) at redshifts z = 9 (left; 
global mass-weighted ionized fraction Xm = 0.28), z = 7 (middle; Xm = 0.94) and z = 6.0 (right; Xm = 0.9999). For reference, the green, dotted line 
shows the result if each source is assumed 50% absorbed, which would be the case if e.g. all of the blue wing of the emission line were absorbed, while all of 
the red wing were transmitted. The error bar in each bin reflects the number of sources in that bin found in our computational volume, (top panels) Bin-by-bin 
ratio of the observed to the intrinsic luminosity function. 



of r_i2 from 10 (for cells close to sources) down to essentially 
(in neutral and shielded cells). There are three broad peaks of the 
distribution, at r_i2 ~ 1 (the H II regions), r_i2 ~ 10~^ (cells at 
and around I-fronts), and r_i2 ~ (neutral regions). By z = 6.0, 
which is well after overlap both the neutral and self- shielded re- 
gions and the I-fronts have mostly disappeared and r_i2 > 0.1 
almost everywhere, rising with time as more galaxies form. These 
values are fairly high compared to the photoionization rate values 
found from the Lyman- a forest at z ^ 2 — 4 (Cen & McDonald 
2OO2I: iTvtler & et alJl2004l : iFan & et alj[200^ : iBolton & Haehnelt 
2007b'), in agreement with the high mean transmitted flux we found 
in ^ 14.31 Both point to somewhat lower ionizing source efficiencies 
than the ones assumed here and to a correspondingly later end of 
reionization. 



4.5 Luminosity function of high-z Ly-a sources 

The luminosity function is an important statistical measure of the 
properties of high=redshift galaxies. It depends on both the intrin- 
sic luminosity of the galaxies and the absorption in the surrounding 
IGM. In Figure [20l we show our results of the size of the effect 
of absorption on a luminosity function of high z objects. For this 
fiducial case we assume that the Ly-a luminosity is simply propor- 
tional to the mass of our halos (similar to our model for the ionizing 
sources). We compute the reduction in luminosity of each halo due 
to absorption (Figures [MTTI show examples of this suppression). 
We assume that the intrinsic Lyman-a emission line is a Gaussian 
with an rms of 160kms~^. Luminosity function of high-redshift 
sources without (in this fiducial case this is just the halo mass func- 
tion; solid line) and with absorption included (dashed) at redshifts 
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e.g. all of the blue wing of the emission line were absorbed, while 
all of the red wing were transmitted. Top panels show the bin-by- 
bin ratios of the observed to the intrinsic luminosity function. Note 
that due to the binning at fixed luminosity (intrinsic or observed) 
this ratio is not the same as the average suppression per source of 
a given mass (i.e. the absorption shifts the curve both down and to 
the left). 

The luminosity function exhibits clear evolution from high to 
low redshift. As we discussed in ^ 14.11 at high redshift the damping 
wings are strong and thus not only the blue side, but also signif- 
icant part of the red side of the line is absorbed, as evidenced by 
the significant difference between the dashed and dotted lines. As a 
result of this absorption, the number of sources per luminosity bin 
drops by one to two orders of magnitude. During the later stages of 
reionization the damping wing effectively disappears. As a conse- 
quence, the change in the faint end of the luminosity function due 
to IGM absorption is on average well-represented by simply reduc- 
ing each source luminosity by 50%, which would be the case if the 
blue half of the line were absorbed and the red half were not. The 
situation is different at the bright end of the luminosity function, 
however, where on average significantly more than half of the in- 
trinsic flux is absorbed at both z = 7 and z = 6. The shape of the 
luminosity function shows some evolution, as well, which is in part 
due to an evolution in the shape of the halo mass function and in 
part to the higher mean absorption for the more massive sources. 
The higher average absorption levels for the luminous sources is 
consequence of the infall which surrounds the high density peaks 
they are in. In order to demonstrate this, we re-calculated the lumi- 
nosity function at 2; = 8.1 and z = 6 using exactly the same data, 
but setting all peculiar velocities to zero. The results are shown in 
Figure|2T| With no peculiar velocities present the intrinsic emission 
of all sources is absorbed on average at roughly the same level, by 
factor of~ lOatz = 8.1 and by factor of 2 at z = 6. The resulting 
luminosity function at late times agrees well with the one where we 
simply assumed 50% absorption. Early-on this is not the case as a 
consequence of the still-present damping wing. 

How important is our assumption that all sources have same 
rms width of their intrinsic emission line profile? To check this, we 
replaced this assumption with one where the rms lin e width varies 
with the halo mass as 133 km/s(M/lO^^M0)^/^ toijkstra et a ll 
l2007h . Results, again at 2; = 6 and in head-to-head comparison 
with our fiducial case of fixed line width, are shown in Figure |22] 
The faint end of the luminosity function proves insensitive to the 
line width, which is easy to understand. As we have shown above, 
on average the IGM completely absorbs the blue half of the line for 
the weaker sources and completely transmits the red half. However, 
the variable line width has some effect on the absorption of bright 
sources. Their lines become wider under this assumption and thus 
are less affected by the absorption due to the infalling gas, resulting 
in higher transmission by a factor of 2. 

The probability distribution of the transmission fraction per 
source of a given mass/luminosity is shown in Fig. |23] There are 
10 random LOS per source and we required at least 50 LOS in 
each bin for sampling the distribution properly. We also plot the 
bin-by-bin average transmission. For the luminosity bins which do 
not contain our minimum number of LOS we plot only the mean. 
Several interesting trends emerge. The distributions are fairly wide 
at all times and for all sources, reflecting the large variations in 
opacity from source to source and from LOS to LOS. The former 



the mean and the median curves are very similar for all bins, and 
are essentially identical at the bright end. 

The distribution itself changes its character as the evolution 
progresses. At early times the distribution is much wider for the 
fainter sources and the mean and median are gently rising towards 
the bright end, reflecting the fact that bright sources are found in 
the middle of larger H II regions, while fainter sources are found 
in a variety of environments. Thus, during the early evolution the 
main factor shaping the distribution is the the local variation of the 
neutral fraction around each source. At late times (z < 7), however, 
the situation changes to the opposite, with the distribution becom- 
ing wider at the bright end and the mean and median decreasing 
there as well. By that time the IGM is already largely ionized and 
the main environmental dependence is due to the anisotropics of 
the density field and, even more so, of the infall around the bright 
peaks, as discussed above. This is clearly demonstrated in Fig. |24l 
where we show the distributions as they would be if there were no 
peculiar velocities. At early times the results are largely unchanged, 
while at late times the variations between the different LOS essen- 
tially disappear and all the curves become flat, showing that the 
distribution is shaped mainly by the effects of the peculiar veloci- 
ties. 

Finally, in Fig.|25]we show the effect of varying intrinsic line 
width on the distributions. Compared to our fiducial case of con- 
stant line width, at high redshift the distributions are hardly af- 
fected, except for slightly higher absorption of the faintest sources. 
However, at later times the varying line width has a more significant 
effects. At 2; = 7 the mean and the median values for the majority 
of sources decrease from ^ 40 — 45% down to ^ 30%, but the 
brightest sources are affected much less. As a result the curves for 
the mean and median become largely flat rather than decreasing to- 
wards the bright end. Furthermore, the distribution becomes wider 
at the faint end, with many more LOS being absorbed by factor of 
10 or more. A similar effect is seen at z = 6, except in this case the 
brightest sources are even less absorbed due to their wider emis- 
sion lines, in agreement with what we observed in the luminosity 
functions. 

Our derived luminosity functions assume a constant mass-to- 
light ratio and are in arbitrary units (halo mass x absorption by the 
IGM), proportional to a yet undetermined mass-to-observed light 
ratio. We can roughly determine the latter by comparing the num- 
ber densities of observed and simulated obiects. lKashikawa & et al.l 
(2006) currently provide the best set of data at 2; > 6. They have 
provided fits to a Schechter function: 

= (^)\.p(-^)^ (2) 

Since the high-redshift data still has large uncertainties, partic- 
ularly in terms of the faint-end slope, the data is fit by assum- 
ing a = (— 2, — 1.5, — 1), with best fit parameters at 2; = 
6.56 given by /op(L* /hf^ferg/s) = (42.74,42.60,42.48) and 
/o^((/)*/Mpc-^h?o) = (-3.14, -2.88, -2.74), respectively. We 
plot these fits in Figure [26] against our derived luminosity function. 
The latter was obtained by rescaling our arbitrary luminosity units 
to physical ones using a constant ratio, 

L = L{Mq) X 10^°-^ (3) 

so as to match it to the observed luminosity function for the same 
number densities of objects. 

The fit assuming a = —1.5 provides by far the best match 
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Figure 21. Luminosity function at redshifts z=8.1 (left; Xm = 0.62) and z=6 (right; Xm = 0.9999) if peculiar velocities are ignored (blue, long-dashed). For 
reference, we also show the data from Fig.[20l same notation. 
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Figure 22. Luminosity function at redshifts z=8.1 (left; Xm = 0.62) and z=6 (right; Xm = 0.9999) for variable emission line width (blue, long-dashed; see 
text for details). For reference, we also show the data from Fig.|20] same notation. 



luminous and the faint end should both be expected, the former due 
to cosmic variance and the latter due to both numerical resolution 
and lack of reliable observational data. Matching the other two faint 
end slopes would require us to relax our constant mass-to-light ratio 
assumption. 

Taking equation [3] at face value, we can now make an approx- 
imate correspondence between observed luminosities and masses 
of the un derlying halos. The sourc es observed with Subaru at 
z = 6.56 lKashikawa & et all J2006h have luminosities IQ^^ _ 



10 Mq or larger, or relatively rare halos. The observations are not 
yet sufficiently sensitive to detect the faint end, which contributes 
most of the ionizing emissivity during reionization. Hence, claims 
that observations show that there are not enough ionizing photons 
at 2 ~ 6 to reionize the universe appear premature. 



5 CORRELATION FUNCTIONS 
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Figure 23. Transmission fraction as a function of luminosity at redshifts z = 9.0 (xm = 0.28; left), z = 7.0 (xm = 0.94; center) and z = 6.0 
(xm = 0.9999; right). Shown are (dashed lines, top to bottom) 0.023, 0.16, 0.5, 0.84, and 0.977 percentiles (e.g., 2.3% of points have suppression less than 
the uppermost line). The center line is the median of all the LOS. We required at least 50 LOS in each bin for sampling the distribution properly. There are 10 
random LOS per source. We also plot the mean (solid line) for all LOS in each bin. 






Figure 24. Same as Fig.|23] but ignoring any peculiar velocities of the halos and the IGM. 






Figure 25. Same as Fig.|23] but for a variable Ly-a emission line width, as discussed in the text. 
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Figure 27. Projection of the sources as seen by a mock flux-limited survey with L > lO^-^M© (198 sources in total) at z = 9 (left panel; Xn 



0.28) 



and sources with the same number density if the IGM absorption were ignored (middle panel) and the 2-point 3D correlation functions (right panels) of the 
distribution with IGM absorption (solid) and without (dashed) and their ratio (top). 



100 




z = 7 



40 60 
Y [Mpc/h] 



40 60 
Y [Mpc/h] 



100 




1 10 
R [cMpc/h] 



Figure 28. Projection of the sources as seen by a mock flux-limited survey with L > 10^^ Mq (1617 sources in total) at z = 7 (left panel; Xn 



0.94) 



and sources with the same number density if the IGM absorption were ignored (middle panel) and the 2-point 3D correlation functions (right panels) of the 
distribution with IGM absorption (solid) and without (dashed) and their ratio (top). 



The latter clustering ha s been recently observed at z ^ 5.7 
dlVluravama & et al]|2007h . An interesting question to ask is if the 
absorption due to the surrounding IGM affects this clustering. If 
this were the case, then measuring the correlation function of high- 
redshift Ly-a sources can give us information about the state of the 
IGM at that time. 

We derive the correlation functions as follows. First we cal- 
culate the total luminosity of each source with and without IGM 
absorption using the same method as above, but instead of ran- 
dom LOS directions we only consider parallel LOS, as would be 
seen by far away observer. For simplicity we consider LOS paral- 
lel to the axes of our computational box. We mock a flux-limited 
survey by imposing a cutoff on the observed luminosity. We com- 
pare the resulting correlation function to the one obtained for the 
same number, but now of the brightest sources based on their in- 
trinsic luminosity (i.e. the ones hosted by the most massive halos, 
thus ignoring IGM absorption in this latter case). We calculate the 
3D correlation func tions by direct summation over all pairs of ha- 
los as described in iMartell (Il99lh . The results at redshift z = 9 
(with cutoff Lmin(M) = lO^-^Mo) and z = 7 (with cutoff 



ence in the correlation functions with and without IGM absorp- 
tion, while at the same time allowing for sufficient number of halos 
above the cutoff to reduce the noise of the correlation. 

In Figure [27] we show the projection at z = 9 of the two 
source distributions onto the Y-Z plane (left), and the correspond- 
ing 3D correlation functions and their ratio (right). The projection 
shows that the two source populations differ, but cluster in the same 
spatial regions. This visual impression is confirmed by the correla- 
tion functions. We find that IGM absorption introduces only small 
variations in the correlation of sources. The difference is largest at 
small scales, for separations below 1 comoving Mpc. Even there 
they never exceed 10%. At intermediate scales, R ^ 2 — 10 Mpc 
the departures are up to 5%. There are no appreciable differences 
at large scales. 

At redshift z = 7, close to overlap (Figure [28]) there are 
many more sources, even with the low luminosity cutoff raised to 
10^° M©. The two source distributions remain different, but are still 
clustered in a very similar way, resulting in largely identical corre- 
lation functions, which never differ by more than 0.7%. 

The small effect of IGM absorption on the correlation function 
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Figure 26. Simulated luminosity f unction at z = 6.6 (dashed; Xm = 0.99) 
vs. best fits of Kashi kawa & et al.l 12006) (at z = 6.56) for (top to bottom) 
faint-end slopes ofa = — 2,-1.5,— 1) (solid). 

strongly clustered and the ionization field is closely correlated with 
the galaxy field. The luminous sources are typically found in the 
inner parts of the largest ionized bubbles, where they are typically 
unaffected by damping from the remaining (few) neutral patches. 
This is also supported by the profiles in Fig. [12] which show the 
damping wing effect being important over the same time interval 
for both luminous and average sources. The damping affects aver- 
age sources more strongly, since these are more often found closer 
to a neutral patch than the more massive sources. However, this 
does not change the correlation function significantly. Basically, 
the reionization patchiness has little effect on the source clustering 
properties. The reason is that while the IGM absorption diminishes 
the flux from all sources, the same source clusters (although not 
necessarily the same individual sources) are seen as would be with- 
out IGM absorption. 

At late times any clearly defined ionized regions have already 
disappeared and the effect of IGM absorption is to replace some 
sources above the luminosity cutoff with other ones essentially at 
random, due to small local variations of the residual neutral fraction 
and the gas velocities. Therefore, the reionization patchiness ulti- 
mately has little effect on the correlation function. That of course 
does not mean that the source population is not modified by the 
IGM - a flux-limited survey will see many fewer sources than if 
IGM absorption were not present, but the clustering properties of 
those sources are almost the same as with no absorption for the 
same number d ensity of source s. 

Recently, McOui nn et claimed that reionization 

patchiness has a significant effect on observed s ource clustering, 
in app a rent contradiction to our resul ts (see also iFurlanetto et aP 
l2006bl : iMesinger & Fu rlanettd l2008b'). However, they compared 
the clustering properties of Ly-a sources with and without (i.e. in- 
trinsic) IGM absorption for a fixed luminosity cutoff, rather than 
at fixed number density of sources, as we did. When the IGM ab- 
sorption is accounted for, the sources remaining above the imposed 



halos. This naturally results in higher bias of the observed sources 
compared to all sources with intrinsic luminosity above that cut- 
off. However, it is not necessarily related to reionization patchi- 
ness. E.g. source bias will increase also if the small residual neutral 
fraction in the ionized IGM increases, boosting the IGM opacity 
and causing dimmer, less clustered sources to fall below the lumi- 
nosity cutoff, thus no neutral patches are required for this to hap- 
pen. In fact, our simulation results show almost no neutral patches 
below z ~ 6.6 (at which time the global mean ionized fraction 
by mass is belo w 1%). Therefore, whil e our conclusions disagree 
with the ones oflM cOuinn et "Zl (l2007h . at least some of the dif- 
ferences could be attributed to the different comparisons we make, 
but so me of the variations are possibly real, since Mc Ouinn et al.l 
(l2007h claimed that the trend is present even for fixed number den- 
sities, albeit with no quantitative details. We would like to stress, 
however, that we observe the same qualitative trend of enhanced 
clustering of hy-a sources due to patchiness during the early stages 
of reionization, but the quantitative level of the effect is different, 
being much weaker in our case. 

There are a number of possible explanations. In particular, in 
our simulations the pre-absorption clustering between ionized re- 
gions and sources appears t o be m ore important than in the sim- 
ulations of iMcOuinn et al ] toot . There are also notable differ- 
ences between our and their modeling of the Ly-a sources and 
the IGM absorption. McOuinn et al.l (l2007l) (as many other current 
studies do) assume complete resonance absorption of the blue side 
of the Ly-a line, and only study the suppression due to the IGM 
damping wing and do not include velocity effects in their analy- 
sis. They also assume a particular duty cycle for their Ly-a emit- 
ters (that o nly 25% of halo s host emitters), which we do not do in 
this work. McOuinn et al.l (2007) also ran simulations with differ- 
ent minimum source halo mass cutoffs, either significantly higher 
(Mmin = 4 X lO^^Mo) or significantly lower (Mmin lO^M©) 
than the one we have here (Mmin = 2.2 x lO^M©). It is possi- 
ble, therefore, that our results are more relevant than these previous 
works, if the low-mass sources missing in our current simulation 
are in fact strongly suppressed during the late stages of reioniza- 
tion due to Jeans-mass filtering. This can only be resolved conclu- 
sively by detailed simulations which actually follow the compli- 
cated radiative feedback effects on low-mass halos, which is well 
beyond the scope of this work. We have run several tests in order 
to try to understand these differences, following suggestions by the 
referee. One possibility we investigated was that our H II regions 
size distribution is more strongly peaked and narrower than the one 
found in the above works, which, if it were the case might have ex- 
plained some of the clustering differences. We found, however, that 
the H II region size distributions deriv ed from our simulation s is in 
fair agreement with the one from the McOuin n et al.l toot sim- 
ulations, and hence this offers no plausible explanation of the dif- 
ferences. We also compared the Ly-a damping wing optical depth 
distributions for sou rces of different mass at a range of reionization 
stages, as derived by IMesinger & Furlanettd (|^08b) (their Figures 
2 and 3) and again found no significant differences between our 
results and theirs. We conclude that more detailed and direct com- 
parisons will be required in order to evaluate and understand any 
differences between our results. 



6 SUMMARY AND CONCLUSIONS 
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tailed structure formation and radiative transfer simulations, which 
allowed us to evaluate many features which can only be studied 
by detailed simulations, as well as to quantify better a number of 
previously-proposed effects. We followed the full reionization his- 
tory self-consistently and accounted for the actual source distribu- 
tion, neutral fraction, density and velocity fields. 

We find that the density, neutral fraction and velocity fields are 
all highly anisotropic, which results in large variations in the IGM 
transmission and source visibility among different LOS. The veloc- 
ity effects, both gas infall and source peculiar velocity are most im- 
portant for massive, luminous sources. The most luminous sources 
are found in highest peaks of the density field, which at late times 
are significantly overdense out to ~ 10 comoving Mpc (cMpc) and 
are surrounded by infall extending to 20 cMpc. The infall of gas 
blueshifts it in frequency space and results in significant absorption 
on the red side of the line center, while the peculiar velocity of the 
source itself can either alleviate or exacerbate this effect, depending 
on the halo- and infall velocity alignment. 

The spherically- averaged local density enhancement and 
gas infall have been modelled analytically in approximate ways 
(IBarkana 2004), and thus can be incorporated in semi-analytical 
models (e.g. iDiikstra et al.ll2007h . However, such models are un- 
able to account for the strong intrinsic anisotropics of the neutral 
fraction, density and velocity fields. The analytical and semianalyt- 
ical models typically assume spherical symmetry and full ioniza- 
tion inside the H II regions, both of which assumptions are quite 
unrealistic. 

The Ly-a lines we derive are generally asymmetric and vary 
hugely from LOS to LOS. The luminous sources form at the high- 
est density peaks and as a consequence their line centers are al- 
ways highly-absorbed even though their proximity regions are very 
highly ionized, with typical neutral fractions xm ^ 10~^ — 10~^. 
The luminous sources also more affected by infall and exhibit more 
pronounced proximity region with higher transmission of the blue 
wing of the line. 

High-redshift sources are strongly clustered around the high 
peaks of the density field. The central source contributes the ma- 
jority of the ionizing flux only in its immediate vicinity, within 1-2 
comoving Mpc. Beyond that distance the ionizing flux is dominated 
by the fainter sources clustered around it. This dominance is partic- 
ularly strong at late times, when both many more sources form and 
the ionized regions become larger, resulting in the fainter sources 
contributing up to 2 orders of magnitude more photons than the 
central source. 

Compared to single- source ionized bubbles, the larger H II re- 
gions from clustered sources diminish the effects from the damping 
wing of the line. Nevertheless, these remain significant until fairly 
late (ionized mass fraction Xm = 0.3 — 0.7, which for the simula- 
tion considered here corresponds to redshifts 2; '-^ 9 — 8). Interest- 
ingly, the average damping wing effect is similar for luminous and 
typical sources, even though naiively one might expect that damp- 
ing could be weaker for the former, since they are typically in the 
middle of large bubbles, away from the neutral patches, unlike the 
fainter sources, which are more evenly distributed. 

Both the mean IGM transmission and the typical photoioniza- 
tion rates we find are high compared to observations at 2; '-^ 6, indi- 
cating that our adopted source efficiencies are also high. The mean 
IGM transmissivity decreases only slowly towards the higher red- 
shifts and the GP transparency occurs significantly after the actual 
overlap epoch. For the simulation considered here overlap (defined 



small fraction (few to 10%) of the flux is transmitted. By overlap 
the spectra start showing significant transmission gaps in the mean 
IGM (i.e. away from the proximity region of a luminous source). 

We find that for a given number density of sources (e.g. as de- 
termined by observations) the clustering of these sources depends 
only weakly on the IGM absorption during reionization. As a con- 
sequence, the reionization patchiness has little effect on the ob- 
served Ly-a source clustering, which implies that source clustering 
is not a good indicator for reionization patchiness. 

Our derived luminosity function assuming constant mass-to- 
light ratio provides an excellent match to the shape of the ob- 
served luminosity function at 2; = 6.6 with faint-end slope of 
a = —1.5. The resulting mass-to-light ratio implies that the ma- 
jority of sources responsible for reionization are too faint to be ob- 
served by the current surveys. 
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